#ifndef __2BWTSEARCH_CPU_CU__
#define __2BWTSEARCH_CPU_CU__

#include"Index.h"
#include "OccValue_cpu.h"
#include "Match_cpu.h"


bool cpu_BackwardSearch(BWT *bwt, unsigned int * Queries, unsigned int queryId, 
		unsigned int lower, unsigned int upper, 
		unsigned int * temp_start, unsigned int * temp_end)
{
	int i;
	char c;
	if(*temp_start > *temp_end)
		return false;

	for (i= (int)(upper-1); i>= (int)lower; i--)
	{
		c = host_getCharFromQueries(Queries, queryId, i);
		if(host_bwtCalc(bwt, c,temp_start, temp_end)== false)
			return false;
	}
	return true;
}




bool cpu_ForwardSearch(BWT* revbwt, unsigned int * Queries, unsigned int queryId, 
		unsigned int lower, unsigned int upper, unsigned int *temp_rev_start, 
		unsigned int *temp_rev_end, unsigned int *temp_start, unsigned int* temp_end)
{
	char c;
	unsigned int occCount[ALPHABET_SIZE];
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];
	unsigned int i;
	int k;

	for (i=lower; i<upper && *temp_start<= *temp_end; i++)
	{
		c = host_getCharFromQueries(Queries, queryId, i);
#if LAYER_INDEX_1
		host_GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, (*temp_rev_start), revbwt->inverseSa0, occCount_start) ;
		host_GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, (*temp_rev_end) + 1, revbwt->inverseSa0, occCount_end) ;
#endif

#if LAYER_INDEX_2
		host_BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, 
						(*temp_rev_start), revbwt->inverseSa0, occCount_start);
		host_BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, 
						(*temp_rev_end)+1, revbwt->inverseSa0, occCount_end);
#endif
		host_bwtCalcWithOccValue(revbwt->countValue, c, occCount_start[c], occCount_end[c], temp_rev_start, temp_rev_end);

		occCount[ALPHABET_SIZE-1] = 0;
		for(k=ALPHABET_SIZE-2; k>=0; k--)
		{
			occCount[k] = occCount[k+1] + occCount_end[k+1] - occCount_start[k+1];
		}

		*temp_end = *temp_end - occCount[c];
		*temp_start = *temp_end - (*temp_rev_end - *temp_rev_start);
	}
	if (*temp_start <= *temp_end && *temp_start > 0)
		return true;
	return false;
}



/**
  @return the number of sa
 */
uint64 cpu_BackwardSearch1Error(uint3* sa,
		BWT * bwt, unsigned int * Queries, unsigned int queryId,
		unsigned int lower, unsigned int upper, unsigned int mismatchCount, 
		unsigned int * sa_start, unsigned int * sa_end, 
		long* numMatch = NULL) 
{
	unsigned int temp_start=(*sa_start),temp_end=*sa_end;
	int i;
	char c, misBase;
	uint64 numSA = 0;


	for (i=upper-1; i >=(int)lower && *sa_start<=*sa_end; i--)
	{
		c = host_getCharFromQueries(Queries, queryId, i);

		for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;
			temp_start = *sa_start;
			temp_end = *sa_end;

			if(!host_bwtCalc(bwt,misBase,&temp_start,&temp_end))
				continue;

			//add pausedQuery to buffer	
			if(cpu_BackwardSearch(bwt,Queries, queryId, 0,i, &temp_start, &temp_end))
			{
				if(temp_end >= temp_start) {
	
						sa[numSA].x = queryId;
						sa[numSA].y = temp_start;
						sa[numSA].z = temp_end;
						numSA++;

					if(numMatch != NULL) {
						(*numMatch) += (temp_end - temp_start + 1);
						if(*numMatch > MAX_NUM_ALIGN_PER_QUERY) {
							*numMatch = -1;
							return 0;
						}
					}
				}
			}
		}
		host_bwtCalc(bwt, c,sa_start,sa_end);
	}

	return numSA;
}


uint64 cpu_ForwardSearch1Error(uint3* sa, 
		BWT * revbwt, unsigned int * Queries, unsigned int queryId, 
		unsigned int lower, unsigned int upper, unsigned int mismatchCount, 
		unsigned int * sa_rev_start, unsigned int * sa_rev_end, 
		unsigned int * sa_start, unsigned int * sa_end, 
		long* numMatch = NULL)
{
	unsigned int temp_start=*sa_start,temp_end=*sa_end, temp_rev_start=*sa_rev_start, temp_rev_end=*sa_rev_end;
	unsigned int occCount[ALPHABET_SIZE];
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];
	unsigned int i;
	int k;
	char c, misBase;
	uint64 numSA = 0;

	for (i=lower; i<upper && *sa_start<= *sa_end; i++)
	{
#if LAYER_INDEX_1
		host_GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, (*sa_rev_start), revbwt->inverseSa0, occCount_start) ;
		host_GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, (*sa_rev_end) + 1, revbwt->inverseSa0, occCount_end) ;
#endif

#if LAYER_INDEX_2
		host_BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, 
						(*sa_rev_start), revbwt->inverseSa0, occCount_start);
		host_BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, 
						(*sa_rev_end)+1, revbwt->inverseSa0, occCount_end);
#endif

		occCount[ALPHABET_SIZE-1] = 0;
		for (k=ALPHABET_SIZE-2; k>=0; k--)
		{
			occCount[k] = occCount[k+1] + (occCount_end[k+1] - occCount_start[k+1]);
		}

		c = host_getCharFromQueries(Queries, queryId, i);
		for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;

			temp_rev_start = *sa_rev_start;
			temp_rev_end = *sa_rev_end;

			host_bwtCalcWithOccValue(revbwt->countValue, misBase,occCount_start[misBase], occCount_end[misBase], 
					&temp_rev_start, &temp_rev_end);

			temp_end = *sa_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if (temp_start > temp_end)
				continue;


			if(cpu_ForwardSearch(revbwt, Queries, queryId, i+1,upper, &temp_rev_start, &temp_rev_end, &temp_start,&temp_end)) 			{
				//ReportResult(mismatchCount+1, temp_start, temp_end, saInterval);
				if(temp_end >= temp_start) {

						sa[numSA].x = queryId;
						sa[numSA].y = temp_start;
						sa[numSA].z = temp_end;
						numSA++;
			

					if(numMatch != NULL) {
						(*numMatch) += (temp_end - temp_start + 1);
						if(*numMatch > MAX_NUM_ALIGN_PER_QUERY) {
							*numMatch = -1;
							return 0;
						}
					}
				}
			}
		}

		host_bwtCalcWithOccValue(revbwt->countValue, c,occCount_start[c], occCount_end[c], sa_rev_start, sa_rev_end);

		*sa_end = *sa_end - occCount[c];
		*sa_start = *sa_end - (*sa_rev_end - *sa_rev_start);
	}

	return numSA;
}

uint64 cpu_BWTExactMatch(uint3* sa, unsigned int * Queries, unsigned int queryId,
		BWT *bwt, BWT *revbwt, long* numMatch)
{
	unsigned int start, end;
	int i;
	char c;
	unsigned long long packedPattern = 0;
	int queryLen = QUERY_LENGTH;
	unsigned int * query = NULL;
	for (i = 0; i < TABLE_SIZE ; i++) {
		packedPattern<<=2;
		c = host_getCharFromQueries(Queries,queryId, queryLen - TABLE_SIZE + i);
		packedPattern |= (c & 3);
	}
	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];

	if(cpu_BackwardSearch(bwt,Queries, queryId, 0, QUERY_LENGTH-TABLE_SIZE, &start, &end)){
		sa->x = queryId;
		sa->y = start;
		sa->z = end;
		if(numMatch != NULL) {
			(*numMatch) += (end - start + 1);
			if(*numMatch > MAX_NUM_ALIGN_PER_QUERY) {
				*numMatch = -1;
				return 0;
			}
		}		
		return 1;
	}
	return 0;
}

uint64 cpu_BWT1MismatchSearch(uint3* sa, unsigned int * Queries, unsigned int queryId,
		BWT *bwt, BWT *revbwt, long* numMatch){

	//int sizeX = (QUERY_LENGTH / 2);
	//int sizeY = QUERY_LENGTH;
	int i;
	unsigned int start,end,revStart,revEnd;
	uint64 numSA = 0;
	unsigned long long l_packedPattern = 0;
	unsigned long long r_packedPattern = 0;
	unsigned long long rev_packedPattern = 0;
	unsigned long long packedPattern = 0;

	for (i = 0; i <bwt->lookupTable->tableSize ; i++) {
		packedPattern<<=2;
		packedPattern |= host_getCharFromQueries(Queries, queryId, QUERY_LENGTH - bwt->lookupTable->tableSize+i) & 3;
	}
	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];
	cpu_BackwardSearch(bwt, Queries, queryId , (QUERY_LENGTH / 2), QUERY_LENGTH-bwt->lookupTable->tableSize, &start, &end);

	if(start <= end)
		numSA += cpu_BackwardSearch1Error(sa + numSA,
				bwt, Queries, queryId,0,(QUERY_LENGTH/2),0,
				&start, &end, numMatch);

	if(numMatch != NULL && (*numMatch) == -1)
		return 0;
	for (i = 0; i < bwt->lookupTable->tableSize ; i++) {
		l_packedPattern<<=2;
		l_packedPattern |= (host_getCharFromQueries(Queries, queryId, i)&3 );
	}
	r_packedPattern = l_packedPattern;
	l_packedPattern <<= (bwt->lookupTable->tableSize-revbwt->lookupTable->tableSize)*2;
	r_packedPattern <<= (bwt->lookupTable->tableSize-revbwt->lookupTable->tableSize)*2;
	r_packedPattern |= (1<<(bwt->lookupTable->tableSize-revbwt->lookupTable->tableSize)*2) - 1;

	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		rev_packedPattern |= (host_getCharFromQueries(Queries, queryId, TABLE_SIZE-i-1) & 3);
	}

	start = l_packedPattern ? bwt->lookupTable->table[l_packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[r_packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	if (end - start < revEnd - revStart) {
		revStart =revEnd - (end - start );
	} else if (end - start > revEnd - revStart) {
		start = end - (revEnd - revStart);
	}

	cpu_ForwardSearch(revbwt, Queries, queryId, revbwt->lookupTable->tableSize, (QUERY_LENGTH/2), &revStart, &revEnd, &start, &end);
	if(start <= end)
		numSA += cpu_ForwardSearch1Error(sa + numSA, revbwt, Queries, queryId, (QUERY_LENGTH/2), QUERY_LENGTH, 0, &revStart, &revEnd, &start, &end, numMatch);
	return numSA;
}
/**
  @return the number of SA
 */
uint64 cpu_BWT2MismatchCaseA(uint3* sa,	
		unsigned int * Queries, unsigned int queryId,  
		BWT *bwt, BWT *revbwt, long* numMatch) {
	unsigned int temp_start=0,temp_end=0;
	unsigned long long packedPattern = 0;
	unsigned int i;
	char c, misBase;
	unsigned int start, end;
	uint64 numSA = 0;


	for (i = 0; i < TABLE_SIZE ; i++) {
		packedPattern<<=2;
		c = host_getCharFromQueries(Queries, queryId, QUERY_LENGTH- TABLE_SIZE +i);
		packedPattern |= (c & 3);
	}
	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];


	if(start <= end &&  cpu_BackwardSearch(bwt, Queries, queryId,
				sizeY, QUERY_LENGTH - TABLE_SIZE, &start, &end))
	{
		for (i=sizeY-1; i>=1 && start<=end; i--)
		{
			c = host_getCharFromQueries(Queries, queryId, i);
			for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
			{
				if( misBase == c)
					continue;
				temp_start = start;
				temp_end = end;
				if(host_bwtCalc(bwt,misBase,&temp_start,&temp_end)){
					numSA += cpu_BackwardSearch1Error(sa + numSA,
							bwt, Queries, queryId,0,i,1, 
							&temp_start, &temp_end, numMatch);

					if((numMatch != NULL) && (*numMatch == -1)) {
						return 0;
					}
				}
			}
			host_bwtCalc(bwt, c, &start, &end);
		}
	}

	return numSA;
}


uint64 cpu_BWT2MismatchCaseB(uint3* sa,	
		unsigned int * Queries, unsigned int queryId,  
		BWT *bwt, BWT *revbwt, long* numMatch) {
	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;
	unsigned int occCount[ALPHABET_SIZE];

	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];

	unsigned long long l_packedPattern = 0;
	unsigned long long r_packedPattern = 0;
	unsigned long long rev_packedPattern = 0;

	//unsigned long long mask;
	int k;
	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;
	int numSA = 0;

	l_packedPattern = 0;
	r_packedPattern = 0;
	rev_packedPattern = 0;
	// look-up the first characters
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		l_packedPattern<<=2;
		c = host_getCharFromQueries(Queries, queryId, i);
		l_packedPattern |= (c & 3);
	}
	r_packedPattern = l_packedPattern;
	// if the look-up tables are of different size
	l_packedPattern <<= (bwt->lookupTable->tableSize -  revbwt->lookupTable->tableSize)*2;
	r_packedPattern <<= (bwt->lookupTable->tableSize - revbwt->lookupTable->tableSize)*2;
	r_packedPattern |= (1<<(bwt->lookupTable->tableSize - revbwt->lookupTable->tableSize)*2) - 1;

	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = host_getCharFromQueries(Queries, queryId, revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3);
	}

	start = l_packedPattern ? bwt->lookupTable->table[l_packedPattern-1]+1 : 1;
	end  = bwt->lookupTable->table[r_packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	//Lookup SA correction...
	if (end - start<revEnd-revStart)
	{
		revStart = revEnd - (end - start);
	} else if (end - start > revEnd - revStart)
	{
		start = end - (revEnd - revStart);
	}

	if (start <= end && cpu_ForwardSearch(revbwt, Queries, queryId, 
				revbwt->lookupTable->tableSize,sizeY,
				&revStart,&revEnd,&start,&end))
	{
		for (i=sizeY; i<sizeZ && start<=end; i++)
		{
#if LAYER_INDEX_1
			host_GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, revStart, revbwt->inverseSa0, occCount_start) ;
			host_GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, revEnd + 1, revbwt->inverseSa0, occCount_end) ;
#endif
#if LAYER_INDEX_2
		host_BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, 
						revStart, revbwt->inverseSa0, occCount_start);
		host_BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, 
						revEnd+1, revbwt->inverseSa0, occCount_end);
#endif
			occCount[3]=0;
			for (k=2;k>=0;k--) {
				occCount[k]=occCount[k+1]+occCount_end[k+1]-occCount_start[k+1];
			}

			c = host_getCharFromQueries(Queries, queryId, i);
			for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
			{
				if(misBase == c)
					continue;
				temp_rev_start = revStart;
				temp_rev_end = revEnd;

				host_bwtCalcWithOccValue(revbwt->countValue,misBase,occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

				temp_end = end - occCount[misBase];
				temp_start = temp_end - (temp_rev_end - temp_rev_start);

				if(temp_start <= temp_end){
					numSA += cpu_ForwardSearch1Error(sa + numSA, revbwt, 
							Queries, queryId,i+1,sizeZ, 1, &temp_rev_start, &temp_rev_end,
							&temp_start, &temp_end, numMatch);

					if((numMatch != NULL) && (*numMatch == -1)) {
						return 0;
					}							
				}
			}

			host_bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

			end = end - occCount[c];
			start = end - (revEnd - revStart);
		}
	}

	return numSA;							
}



uint64 cpu_BWT2MismatchCaseC(uint3* sa,	
		unsigned int * Queries, unsigned int queryId,  
		BWT *bwt, BWT *revbwt, long* numMatch) {
	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;
	unsigned int occCount[ALPHABET_SIZE];

	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];

	unsigned long long packedPattern = 0;
	unsigned long long rev_packedPattern = 0;

	unsigned long long ALLONE = (1<<(bwt->lookupTable->tableSize*2))-1;
	unsigned long long mask;

	int k;
	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;

	packedPattern = 0;
	rev_packedPattern = 0;

	uint64 numSA = 0;

	for (i = 0; i <bwt->lookupTable->tableSize ; i++) {
		packedPattern<<=2;
		c = host_getCharFromQueries(Queries, queryId, i);
		packedPattern |= (c & 3);
	}
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = host_getCharFromQueries(Queries, queryId, revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3 );
	}

	//for error happens in lookup range
	for(i=sizeX; i<bwt->lookupTable->tableSize; i++)
	{
		c = host_getCharFromQueries(Queries, queryId, i);
		for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
		{
			if(misBase == c)
				continue;
			unsigned long long err_packedPattern = packedPattern;
			unsigned long long err_rev_packedPattern = rev_packedPattern;

			unsigned int bitPos = (bwt->lookupTable->tableSize-i-1)*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_packedPattern |= (3 << bitPos);
			err_packedPattern &= mask;

			bitPos = i*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_rev_packedPattern |= (3 << bitPos);
			err_rev_packedPattern &= mask;

			start = err_packedPattern ? bwt->lookupTable->table[err_packedPattern-1]+1 : 1;
			end = bwt->lookupTable->table[err_packedPattern];
			revStart = err_rev_packedPattern ? revbwt->lookupTable->table[err_rev_packedPattern-1]+1 : 1;
			revEnd = revbwt->lookupTable->table[err_rev_packedPattern];

			//Lookup SA correction...
			if (end - start < revEnd - revStart) {
				revStart = revEnd - ( end - start);
			}else if ( end - start > revEnd - revStart) {
				revStart = revEnd -( revEnd - revStart);
			}

			if(start > end)
				continue;
			unsigned int pos = revbwt->lookupTable->tableSize;
			cpu_ForwardSearch(revbwt,Queries, queryId, pos, sizeY, &revStart, &revEnd, &start, &end);

			if( start <= end ){
				numSA += cpu_ForwardSearch1Error(sa + numSA, 
						revbwt,  Queries, queryId,sizeY, sizeZ, 1, 
						&revStart,&revEnd,&start, &end, numMatch);

				if((numMatch != NULL) && (*numMatch == -1)) {
					return 0;
				}	
			}
		}
	}

	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];
	revStart = packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	//Lookup SA correction...
	if (end - start < revEnd - revStart) {
		revStart = revEnd - ( end - start);
	} else if (end - start > revEnd - revStart) {
		start = end - ( revEnd - revStart);
	}

	if(sizeX > revbwt->lookupTable->tableSize){
		cpu_ForwardSearch(revbwt, Queries, queryId, TABLE_SIZE ,sizeX,
				&revStart,&revEnd,&start,&end);
		i = sizeX;
	}
	else 
		i = TABLE_SIZE;
	for (; i<sizeY && start<=end ; i++)
	{
		c = host_getCharFromQueries(Queries, queryId, i);
#if LAYER_INDEX_1	
		host_GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, revStart, revbwt->inverseSa0, occCount_start) ;
		host_GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, revEnd + 1, revbwt->inverseSa0, occCount_end) ;
#endif
#if LAYER_INDEX_2
		host_BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, 
						revStart, revbwt->inverseSa0, occCount_start);
		host_BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, 
						revEnd+1, revbwt->inverseSa0, occCount_end);
#endif

		occCount[ALPHABET_SIZE-1] = 0;
		for(k=ALPHABET_SIZE-2; k>=0; k--)
		{
			occCount[k] = occCount[k+1] + occCount_end[k+1] - occCount_start[k+1];
		}

		for(misBase=0; misBase < ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;
			temp_start = start;
			temp_end = end;
			temp_rev_start = revStart;
			temp_rev_end = revEnd;

			host_bwtCalcWithOccValue(revbwt->countValue,misBase,occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

			temp_end = temp_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if(temp_start > temp_end)
				continue;
			if(cpu_ForwardSearch(revbwt, Queries, queryId, i+1,sizeY,&temp_rev_start,&temp_rev_end,&temp_start,&temp_end) && temp_start <= temp_end) {
				numSA += cpu_ForwardSearch1Error(sa + numSA, 
						revbwt, Queries, queryId,sizeY,sizeZ,1,&temp_rev_start,
						&temp_rev_end,&temp_start, &temp_end, numMatch);
				if((numMatch != NULL) && (*numMatch == -1)) {
					return 0;
				}	
			}
		}

		host_bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

		end = end - occCount[c];
		start = end - (revEnd - revStart);
	}

	return numSA;
}



uint64 cpu_BWT2MismatchCaseD(uint3* sa,	
		unsigned int * Queries, unsigned int queryId,  
		BWT *bwt, BWT *revbwt, long* numMatch) {

	uint64 numSA = 0;

	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;
	unsigned int occCount[ALPHABET_SIZE];
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];
	unsigned long long packedPattern = 0;
	unsigned long long rev_packedPattern = 0;

	unsigned long long ALLONE = (1<<(TABLE_SIZE * 2))-1;
	unsigned long long mask;

	int k;
	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;

	packedPattern = 0;
	rev_packedPattern = 0;

	for (i = 0; i < TABLE_SIZE ; i++) {
		packedPattern<<=2;
		c = host_getCharFromQueries(Queries, queryId, sizeX+i);
		packedPattern |= (c & 3);
	}
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = host_getCharFromQueries(Queries, queryId, sizeX+revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3);
	}

	//For error happen in lookup range....
	for (i = sizeY - sizeX; i< TABLE_SIZE ;i++) {
		c = host_getCharFromQueries(Queries, queryId, sizeX + i);
		for (misBase=0;misBase<4;misBase++) {
			if (misBase == c) continue;
			unsigned long long err_packedPattern = packedPattern;
			unsigned long long err_rev_packedPattern = rev_packedPattern;
			unsigned int bitPos = (TABLE_SIZE - i-1)*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_packedPattern |= (3 << bitPos);
			err_packedPattern &= mask;

			bitPos = i*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_rev_packedPattern |= (3 << bitPos);
			err_rev_packedPattern &= mask;

			start = err_packedPattern ? bwt->lookupTable->table[err_packedPattern-1]+1 : 1;
			end = bwt->lookupTable->table[err_packedPattern];
			revStart = err_rev_packedPattern ? revbwt->lookupTable->table[err_rev_packedPattern-1]+1 : 1;
			revEnd = revbwt->lookupTable->table[err_rev_packedPattern];

			//Lookup SA correction...
			if (end - start < revEnd - revStart) {
				revStart = revEnd - ( end - start );
			} else if ( end - start >revEnd - revStart) {
				start = end - (revEnd - revStart);
			}
			if( start <= end && cpu_ForwardSearch(revbwt, Queries, queryId, sizeX+revbwt->lookupTable->tableSize, QUERY_LENGTH, 
						&revStart, &revEnd, &start, &end )) {
				numSA += cpu_BackwardSearch1Error(sa + numSA, bwt, Queries, queryId, 
						0, sizeX, 1, &start,&end, numMatch);
				if((numMatch != NULL) && (*numMatch == -1)) {
					return 0;
				}	
			}
		}
	}

	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	//Lookup SA correction...
	if (end-start<revEnd-revStart) {
		revStart=revEnd-(end-start);
	} else if (end-start>revEnd-revStart) {
		start=end-(revEnd-revStart);
	}

	if(sizeY - sizeX > TABLE_SIZE) {
		cpu_ForwardSearch(revbwt, Queries, queryId, sizeX + TABLE_SIZE ,sizeY,
				&revStart,&revEnd,&start,&end);
		i = sizeY;

	}
	else
		i = sizeX + TABLE_SIZE;
	for (; i<sizeZ && start<=end; i++)
	{
		c = host_getCharFromQueries(Queries, queryId, i);
#if LAYER_INDEX_1
		host_GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, revStart, revbwt->inverseSa0, occCount_start) ;
		host_GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, revEnd + 1, revbwt->inverseSa0, occCount_end) ;
#endif

#if LAYER_INDEX_2
		host_BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, 
						revStart, revbwt->inverseSa0, occCount_start);
		host_BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, 
						revEnd+1, revbwt->inverseSa0, occCount_end);
#endif
		occCount[ALPHABET_SIZE-1] = 0;
		for(k=ALPHABET_SIZE-2; k>=0; k--)
		{
			occCount[k] = occCount[k+1] + occCount_end[k+1] - occCount_start[k+1];
		}

		for (misBase = 0; misBase<ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;
			temp_start = start;
			temp_end = end;
			temp_rev_start = revStart;
			temp_rev_end = revEnd;

			host_bwtCalcWithOccValue(revbwt->countValue,misBase,occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

			temp_end = temp_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if (temp_start <= temp_end)
				cpu_ForwardSearch(revbwt, Queries, queryId, i+1,sizeZ, &temp_rev_start, &temp_rev_end, &temp_start, &temp_end);
			if(temp_start <= temp_end){
				numSA += cpu_BackwardSearch1Error(sa + numSA, bwt, Queries, queryId, 0, 
						sizeX,1,&temp_start, &temp_end, numMatch);
				if((numMatch != NULL) && (*numMatch == -1)) {
					return 0;
				}	
			}

		}
		host_bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

		end = end - occCount[c];
		start = end - (revEnd - revStart);
	}


	return numSA;
}



#endif /*__2BWTSEARCH_CPU_CU__*/
